
// Figure B4

/* =============================================================================
   Step 1: Data processing
   ===========================================================================*/
* Loading BSU link to BSU number in 2014
use "${datain}grk_link_to_grk_number_in_2014", clear 
duplicates drop grk, force
tempfile grklink 
save `grklink', replace

* Loading congestion data for Oslo and Viken from database file 
import dbase using "${datain}losdata_viken" , clear
drop if TID_MRUSH == 0
drop if TID_MRUSH == .
* Averaging morning and afternoon rush
gen TID_RUSH = (TID_MRUSH + TID_ERUSH) / 2

* Measure of congestion: Relative increase in driving time in rush vs outside rush
gen congestion = (TID_RUSH-TID_LAV)/TID_LAV
rename FRA from // Basic statistical unit (grunnkrets) for origin
rename TIL to // Basic statistical unit (grunnkrets) for destination
keep from to congestion
duplicates drop from to, force

* Adding the 2014 BSU numbers to origin BSU
gen grk = from 
merge m:1 grk using `grklink', keep(match master)
drop _merge grk
rename grk2014 from2014 
* Adding the 2014 BSU numbers to destination BSU
gen grk = to 
merge m:1 grk using `grklink', keep(match master)
drop _merge grk
rename grk2014 to2014 

* In case the BSU numbers for 2014 are missing, use existing number
* This is the case for all BSUs that have not changed number over time
replace from = from2014 if from2014 != . 
replace to = to2014 if to2014 != .
drop from2014 to2014
duplicates drop from to, force

* Store dataset
tempfile congestion 
save `congestion', replace

* === Go to main data ==========================================================
use "${dataout}MainDataset" , clear

keep if year >= 2015 & year <= 2017

gen from = grkrets_num 
gen to = grk_bed1_num 

* merging congestion to work for person 1
merge m:1 from to using `congestion', keep(match master)
drop _merge 
rename congestion congestion_to1

rename from temp 
rename to from 
rename temp to 

* merging congestion from work for person 1
merge m:1 from to using `congestion', keep(match master)
drop _merge 
rename congestion congestion_from1

drop from to

gen from = grkrets_num 
gen to = grk_bed2_num 

* merging congestion to work for person 2
merge m:1 from to using `congestion', keep(match master)
drop _merge 
rename congestion congestion_to2

rename from temp 
rename to from 
rename temp to 

* merging congestion from work for person 2
merge m:1 from to using `congestion', keep(match master)
drop _merge 
rename congestion congestion_from2
drop from to

* For each household member ...
forvalues i = 1/2 {
	/* Average of to/from work */
	gen congestion`i' = (congestion_from`i' + congestion_to`i') / 2
	// if congestion_to is missing, use congestion_from as proxy
	replace congestion`i' = congestion_from`i' if congestion`i' == .
	// if congestion_from is missing, use congestion_to as proxy
	replace congestion`i' = congestion_to`i' if congestion`i' == .
}

/* Average value across household members */
egen congestion = rowmean(congestion1 congestion2)

/* Sample selection for regressions */
keep if congestion != . /* Oslo and Viken */
keep if couple == 1 /* only couples */


*============================================================================
// Absorb 3 way FE from variables -> use in binscatters 
*============================================================================

local xvarlist congestion toll_fam_mean ptl_fam_km_mean  

foreach x in `xvarlist' {
	capt drop res_`x'
	capt drop resM_`x'
	* Regressing each variable on geographic fixed effects, storing residuals
	reghdfe `x'  i.year , absorb($FE)  residuals(res_`x') 
	sum `x'  if e(sample)==1
	return list
	dis "`r(mean)'"
	* Creating new variable - residuals + means
	gen resM_`x'=  res_`x' + `r(mean)'
	compress
}

*============================================================================
// HH char: Absorb 3 way FE + HH and ind-work controls from variables  
// Takes a long time to run
*============================================================================

*=========== TOLL ==================

//  congestion   +  toll_fam_mean
local yvarlist    congestion 
local xvarlist    congestion toll_fam_mean
foreach y in `yvarlist' {
	foreach x in `xvarlist' {
		capt drop res_`y'_`x'
		capt drop resM_`y'_`x'
		reghdfe `x'  ptl_fam_km_mean   $age $distance  $time  $publictime  $publicquality   , ///
		absorb($FE   $household  $income  $employment  $education  i.year   ) residuals(res_`y'_`x') 
		sum `x'  if e(sample)==1
		return list
		dis "`r(mean)'"
		gen resM_`y'_`x'=  res_`y'_`x' + `r(mean)'
	}
}


*=========== PTL ==================

// congestion   +  ptl_fam_mean
local yvarlist    congestion 
local xvarlist    congestion ptl_fam_km_mean
foreach y in `yvarlist' {
	foreach x in `xvarlist' {
		capt drop resP_`y'_`x'
		capt drop resPM_`y'_`x'
		reghdfe `x'  toll_fam_mean   $age    $distance  $time  $publictime  $publicquality   , ///
		absorb($FE   $household  $income  $employment  $education  i.year   ) residuals(resP_`y'_`x') 
		sum `x'  if e(sample)==1
		return list
		dis "`r(mean)'"
		gen resPM_`y'_`x'=  resP_`y'_`x' + `r(mean)'
	}
}


* Labeling values (for y-axis on figures)
label var toll_fam_mean "Road toll (NOK)"
label var ptl_fam_km_mean "Bus lane (km)"
label var congestion "Congestion level"
*===============================================================================
// Plot (residualized) correlations between TOLL and HH characteristics
*===============================================================================

*=== loop over the various x charactersitics
local yvarlist toll_fam_mean   		
local xvarlist congestion    
foreach y in `yvarlist' {
	foreach x in `xvarlist' {
		local laby: variable label `y'
		local labx: variable label `x'

		*=== No controls, lfit
		sum `x', det
		capt gen helpvar=1
		reghdfe `y'  `x'  if `x'>=`r(p1)' & `x'<=`r(p99)', ///
			absorb(helpvar)  vce(cluster $clustervar) 
		eststo
		drop helpvar 

		local coeff _b[`x']
		local fmtcoeff : dis %7.5f `coeff'
		dis `fmtcoeff'

		local se _se[`x']
		local fmtse : dis %7.6f `se'
		dis `fmtse'

		sum `x', det
		binscatter  `y'  `x'  if `x'>=`r(p1)' & `x'<=`r(p99)'  , ///
		line(lfit)   nquantiles(20) legend(off)  scale(1.6) lcolor(gs10)   ///
		ytitle("`laby'") ylab(0(10)30)   ///
		xtitle("`labx'")  ///
		graphregion(color(white))   yscale(titlegap(1.5)) xscale(titlegap(1.5)) ///
		subtitle("Coeff: `fmtcoeff'" "(`fmtse')") name(figB4a, replace)
		graph export     "${figures}FigureB4a.png" , replace
		graph export     "${figures}FigureB4a.pdf" , replace
		graph save       "${figures}FigureB4a.gph" , replace

		*=== 3FE, lfit
		sum `x', det
		capt gen helpvar=1
		reghdfe resM_`y'  resM_`x' if `x'>=`r(p1)' & `x'<=`r(p99)' , ///
			absorb(helpvar)  vce(cluster $clustervar) 
		eststo
		drop helpvar 

		local coeff _b[resM_`x']
		local fmtcoeff : dis %7.5f `coeff'
		dis `fmtcoeff'

		local se _se[resM_`x']
		local fmtse : dis %7.6f `se'
		dis `fmtse'

		sum `x', det
		binscatter  resM_`y'  resM_`x'   if `x'>=`r(p1)' & `x'<=`r(p99)' , ///
		linetype(lfit)    nquantiles(20)   legend(off) scale(1.6) lcolor(gs10)  reportreg  ///
		ytitle("`laby'")  ylab(0(10)30)   ///
		xtitle("`labx'")  ///
		graphregion(color(white))   yscale(titlegap(1.5)) xscale(titlegap(1.5)) ///
		subtitle("Coeff: `fmtcoeff'" "(`fmtse')") name(figB4b, replace)
		graph export     "${figures}FigureB4b.png" , replace
		graph export     "${figures}FigureB4b.pdf" , replace
		graph save       "${figures}FigureB4b.gph" , replace
	}
}


*===================================================================================== 
// Plot (residualized) correlations between TOLL and HH characteristics
// Absorb controls using regression! allows for same controls as main regression
*=====================================================================================

*=== loop over the various x charactersitics
local yvarlist    toll_fam_mean    
local xvarlist    congestion    
foreach y in `yvarlist' {
	foreach x in `xvarlist' {
		local laby: variable label `y'
		local labx: variable label `x'

		*=== HH and ind work controls - from detailed regression! , lfit
		sum `x', det
		capt gen helpvar=1
		reghdfe resM_`x'_`y'  resM_`x'_`x' if `x'>=`r(p1)' & `x'<=`r(p99)' ///
			, absorb(helpvar)  vce(cluster $clustervar) 
		eststo
		drop helpvar 

		local coeff _b[resM_`x'_`x']
		local fmtcoeff : dis %7.5f `coeff'
		dis `fmtcoeff'

		local se _se[resM_`x'_`x']
		local fmtse : dis %7.6f `se'
		dis `fmtse'

		sum `x', det
		binscatter  resM_`x'_`y'  resM_`x'_`x'   if `x'>=`r(p1)' & `x'<=`r(p99)' , ///
		linetype(lfit)  nquantiles(20)  legend(off) scale(1.6) lcolor(gs10)  reportreg  ///
		ytitle("`laby'") ylab(0(10)30)  ///
		xtitle("`labx'")  ///
		graphregion(color(white))   yscale(titlegap(1.5)) xscale(titlegap(1.5)) ///
		subtitle("Coeff: `fmtcoeff'" "(`fmtse')")  name(figB4c, replace)
		graph export     "${figures}FigureB4c.png" , replace
		graph export     "${figures}FigureB4c.pdf" , replace
		graph save       "${figures}FigureB4c.gph" , replace
	}
}



*===============================================================================
// Plot (residualized) correlations between PTL and HH characteristics 
*===============================================================================


*=== loop over the various x charactersitics
local yvarlist ptl_fam_km_mean   		
local xvarlist congestion    
foreach y in `yvarlist' {
	foreach x in `xvarlist' {
		local laby: variable label `y'
		local labx: variable label `x'

		*=== No controls, lfit
		sum `x' , det
		capt gen helpvar=1
		reghdfe `y'  `x'  if `x'>=`r(p1)' & `x'<=`r(p99)' ///
			, absorb(helpvar)  vce(cluster $clustervar) 
		eststo
		drop helpvar 

		local coeff _b[`x']
		local fmtcoeff : dis %7.5f `coeff'
		dis `fmtcoeff'

		local se _se[`x']
		local fmtse : dis %7.6f `se'
		dis `fmtse'

		sum `x', det
		binscatter  `y'  `x'  if `x'>=`r(p1)' & `x'<=`r(p99)'  , ///
		line(lfit)   nquantiles(20) legend(off)  scale(1.6) lcolor(gs10)   ///
		ytitle("`laby'") ylab(0(1)3)   ///
		xtitle("`labx'")  ///
		graphregion(color(white))   yscale(titlegap(1.5)) xscale(titlegap(1.5)) ///
		subtitle("Coeff: `fmtcoeff'" "(`fmtse')") name(figB4d, replace)
		graph export     "${figures}FigureB4d.png" , replace
		graph export     "${figures}FigureB4d.pdf" , replace
		graph save       "${figures}FigureB4d.gph" , replace

		*=== 3FE, lfit
		sum `x', det
		capt gen helpvar=1
		reghdfe resM_`y'  resM_`x' if `x'>=`r(p1)' & `x'<=`r(p99)', ///
			absorb(helpvar)  vce(cluster $clustervar) 
		eststo
		drop helpvar 

		local coeff _b[resM_`x']
		local fmtcoeff : dis %7.5f `coeff'
		dis `fmtcoeff'

		local se _se[resM_`x']
		local fmtse : dis %7.6f `se'
		dis `fmtse'

		sum `x', det
		binscatter  resM_`y'  resM_`x'   if `x'>=`r(p1)' & `x'<=`r(p99)' , ///
		linetype(lfit)    nquantiles(20)   legend(off) scale(1.6) lcolor(gs10)  reportreg  ///
		ytitle("`laby'") ylab(0(1)3)   ///
		xtitle("`labx'")  ///
		graphregion(color(white))   yscale(titlegap(1.5)) xscale(titlegap(1.5)) ///
		subtitle("Coeff: `fmtcoeff'" "(`fmtse')") name(figB4e, replace)
		graph export     "${figures}FigureB4e.png" , replace
		graph export     "${figures}FigureB4e.pdf" , replace
		graph save       "${figures}FigureB4e.gph" , replace
	}
}



*===================================================================================== 
// Plot (residualized) correlations between PTL and HH characteristics
// Absorb controls using regression! allows for same controls as main regression
*=====================================================================================


*=== loop over the various x charactersitics
local yvarlist    ptl_fam_km_mean    
local xvarlist    congestion            
foreach y in `yvarlist' {
	foreach x in `xvarlist' {

		local laby: variable label `y'
		local labx: variable label `x'

		* === HH and ind work controls - from detailed regression! , lfit
		sum `x', det
		capt gen helpvar=1
		reghdfe resPM_`x'_`y'  resPM_`x'_`x' if `x'>=`r(p1)' & `x'<=`r(p99)', ///
			absorb(helpvar)  vce(cluster $clustervar) 
		eststo
		drop helpvar 

		local coeff _b[resPM_`x'_`x']
		local fmtcoeff : dis %7.5f `coeff'
		dis `fmtcoeff'

		local se _se[resPM_`x'_`x']
		local fmtse : dis %7.6f `se'
		dis `fmtse'

		sum `x', det
		binscatter  resPM_`x'_`y'  resPM_`x'_`x'   if `x'>=`r(p1)' & `x'<=`r(p99)' , ///
			linetype(lfit)  nquantiles(20)  legend(off) scale(1.6) lcolor(gs10)  reportreg  ///
			ytitle("`laby'") ylab(0(1)3)   ///
			xtitle("`labx'")  ///
			graphregion(color(white))   yscale(titlegap(1.5)) xscale(titlegap(1.5)) ///
			subtitle("Coeff: `fmtcoeff'" "(`fmtse')") name(figB4f, replace)
		graph export     "${figures}FigureB4f.png" , replace
		graph export     "${figures}FigureB4f.pdf" , replace
		graph save       "${figures}FigureB4f.gph" , replace
	}
}


